home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / exp.c < prev    next >
C/C++ Source or Header  |  1990-11-23  |  2KB  |  64 lines
  1. /* EXP and hyperbolic trig functions for Sozobon C.            */
  2. /* Copyright © David Brooks, 1989 All Rights Reserved            */
  3.  
  4. #include <fplib.h>
  5. #include <errno.h>
  6.  
  7. /* EXP:                                    */
  8. /*                                    */
  9. /* Uses the synthesis (2**n) * (2**[m/8]) * b**g            */
  10. /* where b = 2**(1/8) and 8*n+m+g = arg/ln(b).  b**g has a cubic     */
  11. /* approximation: source and accuracy unknown.                */
  12. /*                                    */
  13. /* Beware a bug in the standard release: can't compare two negative    */
  14. /* floating numbers.                            */
  15.  
  16. float exp(a) register float a;
  17. {    fstruct        res;
  18.     register int    aint;
  19.  
  20.     static fstruct huge = {HUGE_AS_INT};
  21.     static float powtab[] = {1.0, 1.09050773, 1.1892071, 1.2968396,
  22.             1.41421356, 1.5422108, 1.68179283, 1.8340081};
  23.  
  24.     if (a > 43.5)
  25.     {    errno = ERANGE;
  26.         return huge.f;
  27.     }
  28.  
  29.     if ((a + 43.5) < 0.0)
  30.         return 0.0;                /* Underflow */
  31.  
  32.     if ((aint = (int)(a *= 11.5415603)) < 0)    /* 8/ln(2) */
  33.         --aint;                /* Correct mathematically */
  34.     a = (float)aint - a;            /* -(frac part)        */
  35.  
  36.     res.f = 1.0 + a * (-0.0866439378 + a * (0.003750577 + 
  37.                 a * -0.11321783e-3));
  38.     res.sc[3] += aint >> 3;            /* Mult by 2**n        */
  39.     return res.f * powtab[aint&7];        /* Mult by 2**[m/8]    */
  40. }
  41.  
  42. /* HYPERBOLIC FUNCTIONS: virtually free                    */
  43.  
  44. float sinh(a) float a;
  45. {    register float    ea = exp(a);
  46.  
  47.     return 0.5 * (ea - (1.0 / ea));
  48. }
  49.  
  50. float cosh(a) float a;
  51. {    register float    ea = exp(a);
  52.  
  53.     return 0.5 * (ea + (1.0 / ea));
  54. }
  55.  
  56. float tanh(a) float a;
  57. {    register float    e2a;
  58.  
  59.     e2a = exp(a);
  60.     e2a = e2a * e2a;            /* exp-squared-a */
  61.     return (e2a - 1.0) / (e2a + 1.0);
  62. }
  63.  
  64.